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Abstract 

Using an improved version of the projection quantum Monte Carlo technique, 
we study the square-lattice Hubbard model with nearest-neighbor hopping t 
and next-nearest-neighbor hopping t' , by simulation of lattices with up to 
20x20 sites. For a given R = 2t'/t, we consider that filling which leads to 
a singular density of states of the noninteracting problem. For repulsive in- 
teractions, we find an itinerant ferromagnet (antiferromagnet) for R = 0.94 
{R = 0.2). This is consistent with the prediction of the T-matrix approxima- 
tion, which sums the most singular set of diagrams. 
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The understanding of itinerant ferromagnetism (FM) is a long-standing problem of solid- 
state physics [0]. In search for a generic model of FM, the Hubbard model, describing 
electrons from a single band subject to a local electron-electron repulsion U, has been inves- 
tigated extensively. Motivated by an exact result of Nagaoka 0, most papers studied the 
stability of a fully polarized state (believed to occur in the phase diagram of the Hubbard 
model at large interaction strength and close to half filling) against single spin flips. The 
results turned out to be strongly dependent on the quasiparticle spectrum and, gener- 
ically, unrealistically large U were required to stabilize a fully polarized state. Motivated 
by the recent proofs of FM in certain models with flat bands |Q, in this Letter we take a 
complementary route, investigating the stability of the paramagnetic phase of a model with 
a high density of states against FM at weak coupling. 

We consider electrons on a square lattice with L = I x I (even /) sites, described by the 
Hamiltonian 

i 

where t and t' are nearest and next-nearest neighbor hoppings, respectively. The bare 
dispersion is = — 2t(cos + cos ky) + At' cos k^ cos ky and its saddle points are, for < 
R = 2t'/t < 1 studied here, at (tt, 0) and (0, tt). In what follows, we always consider the case 
when the chemical potential is — 4t', so that the noninteracting Fermi surface crosses the 
saddle points and the noninteracting density of states at the Fermi level diverges. At this 
so-called Van Hove (VH) density, both the particle-hole and particle-particle susceptibilities 
Xph, Xpp diverge at q = and Q = (vr, tt) and their singular parts are 

Xph(O) = (l/27r2t)(l/v/r^) \n{t/u), 
Xph(Q) = {l/2nh) ln[(l + VT^)/R] ln(t/cu), 

Xpp(O) = {l/4nh){l/Vl^)\n\t/u;), 
Xpp(Q) = {l/2TTH)[aTctan{R/VT^n^)/R]\n{t/iu), 

where u is an infrared energy cutoff. It follows that, in the mean field approximation, 
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infinitesimal interactions cause symmetry breaking: U > implies FM for R > 0.55 and 
antiferromagnetism (AFM) for R < 0.55 and U < leads to superconductivity [^]. The 
magnetic states are interesting because of their metallic character. However, particularly in 
the case of FM, it is known that the use of a bare susceptibility in the Stoner criterion leads 
to incorrect predictions. In order to go beyond mean field theory, let us calculate the U"^ 
correction to the energy of a FM state with magnetization m = {N^ — Ni)/L where N^j is 
the total number of electrons with spin a =1, |, 



with momentum (energy) transfers q ~ 0, Q (cj ~ 0) are cut off in the FM state. We have 
calculated the energy to order f/^ for R = 0.94 on a 64x64 lattice with 1612 electrons [j^ for 
10 closed-shell configurations with magnetizations m < 186/64^ and found that the lowest 
possible m = 2/64^ was stable for all U. Fig. 1 shows a similar result for a 16x16 lattice. 
Thus a resummation of an infinite set of diagrams is necessary. 

In the following, we use the T-matrix approximation (TMA) in order to go beyond 
perturbation theory. TMA is, in general, justified if the interaction range is much shorter 
than interparticle spacing. This situation is realized in the limit R ^ 1. Moreover, scattering 
in the particle-particle channel is the most singular process for all i? > 0. TMA has been 
shown previously to restrict severely the possibility of low-density FM . Following Ref . 
we find that, in TMA, the spin-antisymmetric Landau interaction function /^j^, = —U/[l + 
t/Xpp(k + k', u; = 0)]. This means that the interaction between those points of the Fermi 
surface which are in the neighborhood of the VH points vanishes. Assuming that the shape of 
the noninteracting Fermi surface is not changed, let us consider now, within Landau's Fermi 
liquid theory, the stability of the paramagnetic state against Pomeranchuk instabilities. The 
change of the energy under a local shift Ak o- of the chemical potential for spin-cx electrons 
close to the point k on the Fermi surface is 




m 7^ reduces the always negative E^'^\ h 



ecause the anomalous contributions of processes 
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+ 2 ^ / "^^^^ / c?^'^k'/k:k' Ak,.Ak,,.., (2) 

ct,ct' 

where A^k l/'^k is the angle-resolved density of states, fk is the Fermi velocity, and / dk 
denotes an integral along the Fermi surface. For a ferromagnet Ak,cr = crA and, conse- 
quently, 5E oc A^[l + Fq], where = § dkN]^ § dk' N]^,! f^y / § dkN]^. Let us assume further 
that fk oc 5k, where 5k is the distance between k and the closest VH point. For k, k' both 
in the neighborhood of the same VH point, we expect Xpp(k + k') oc ln^[l/ max(5A;, 
The contribution of such k, k' to the integral for Fq is Fq oc —\ju/t. Thus, intra- 
saddle point processes do not lead to FM at weak coupling. However, if k, k' lie close 
to different VH points, we expect Xpp(k + k') oc ln[l/max(5A;, 5A;')]. This contributes 

Fg {U/t) J^^idk/[k{l + Aln{l/k))] oc -ln[Aln/] where A oc U/t and / is the system 

size. Thus the FM instability criterion Fq < —1 is satisfied for arbitrarily small [/ if / ^ oo. 
Note that TMA predicts a much weaker divergence of F^ with /, as compared to the mean 
field prediction F^ oc —A In/. This is shown explicitly for a 16x16 lattice in Fig. 1. 

The above assumptions that the shape of the Fermi surface does not change and that 
fk oc 5k close to a VH point, are supported by the Quantum Monte Carlo (QMC) data 
for R = 0.94 and f//t = 2 on a 16 x 16 lattice, which shows no appreciable change of the 
Fermi surface. Moreover, our QMC results suggest that the smearing of the momentum 
distribution function n(k) does not increase close to the VH points. This is in contrast to 
the predictions of the second-order perturbation theory in U that the quasiparticle residue 
Z vanishes logarithmically at the VH points, while being finite away from them. 

Remarkably, for electrons described by ^k = k^ky with the Fermi level e = 0, Z is 
finite even at the VH point in TMA. Due to particle-hole symmetry, there is no change 
of the Fermi surface for this model and the imaginary part of the electron self-energy at 
k = is ~ — |£:|/ln^ |£:| for e ^ 0. Thus the wavefunction renormalization Z~^ = 

1 — dTi'/duj\^=Q = 1 — (2/7r) Jl^ deTi" {e) / e"^ does not vanish, although the inverse lifetime 
differs only marginally from the golden-rule result. Anyway, this anomalous lifetime is not 
likely to appear in the transport properties of the VH system. [|l^] 
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The projection QMC method is obtained by filtering out the ground-state component 
from a given trial state = lim.T-^oo^'^'^l'ipT)- Since the imaginary time propagation 
conserves the symmetries of the trial function \iPt), it is also possible to study, e.g., different 
spin subspaces of H by choosing a definite spin for \iPt)- The many-body propagator e~^'^ 
is implemented by use of the one-body time-dependent propagators f/o-(r, 0) defined by the 
discrete Hubbard- Stratonovich fields a{r, r) defined on each site r. We use the usual Trotter 
discretization of the total imaginary time (3 for <T < f3, namely e"^^ = Ea Ua{r, 0). lTl| 
The ground state energy is obtained by Eq = lim Ej3 where 

_ S.(^T|g^.(/^,0)|^T) 

^ j:Ai^T\u^W,o)\iJT) ■ 

Analogously, all relevant correlation functions are obtained by Monte Carlo sampling over 
the fields a{r, r), which become computable provided the trial function \iPt) is chosen to be 
a Slater determinant IS"), as Ucr{T,0)\ipT) remains a Slater determinant too. ^1 
The main improvements of the present scheme are : 

1) Optimization of the trial function \iPt), which has been extended to include not only 
simple Slater determinants IS*) but also the more appropriate Gutzwiller wavefunction {tpg) = 
e~^^\S), with a variational parameter g controlling the total number of doubly occupied sites 
D. This is obtained upon a slight change of the one-body propagator f/<j(r, 0) Ua,ag{j, 0), 
with extra discrete fields (Xgir) needed to decouple the correlated part of the wavefunction 
e^^^ at the initial and final imaginary time r = and r = j3. 

2) Use of the Gutzwiller wavefunction to implement an importance sampling strategy to 
reduce statistical fluctuation of the ground state energy, similar to what is already known 
for the Green function QMC. [1^ We define an estimator that satisfies the "zero variance 
principle", namely that the variance of the energy is zero if the trial wavefunction \iI)t) 
approaches the ground state |V'o)- In order to satisfy this principle the energy has to be 
computed at the initial and final imaginary times (r = 0,/5) since only at those points the 
value of the estimator is independent from the random fields a in the limit when \iPt) iV'o)- 
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We define therefore 

= ' 

wfiere Wo-.o-g = |('S'|f/cr,(Tg(T, 0)|5')| is tlie conventional weiglit (tlie sign being Sg-^o-g)- Tfie 
estimator is 

p (g|g,^...g(AO) + ^.,.g(AO)i/-,|g) 

Ea,.g = ^ , (3) 

and Hg = e^^He~^^ is a nonunitary transformation of the Hamiltonian, which can be 
computed by standard algebra. Suppose now that the ground state is well approximated by a 
Gutzwiller wavefunction (the method can be clearly generalized to any Jastrow wavefunction, 
\'^|Jg)=e~^^^\S) , with J any two-body interaction term), then it is easy to show that the above 
estimator will acquire a small variance, because the left eigenvector of Hg and the right one 
of H_g are well approximated by a single Slater determinant IS"). 

In order to minimize the finite-size effects and also to stabilize the simulation we restrict 
the Slater determinant to a product of plane-wave Slater determinants in each spin sector 
\S) = \Si) ^\Si) , with the condition to fill all degenerate single-particle levels in both spin 
sectors (the closed-shell condition). In this case the total spin is S = \N-^ — Ni\/2. Fig. 2a 
shows that the improvement |]TE[ in the energy estimator allows to obtain accurate values 



of energy in each spin sector even for small imaginary time f3, before the sign problem 
becomes serious. Among the states obtained with a trial state lipx) consistent with the 
closed-shell condition, on a 16x16 lattice the highest spin sector has the lowest energy, as 
shown explicitly by comparing to the lowest spin state (singlet in the thermodynamic sense). 
We emphasize that, for a finite lattice, the instability occurs at a finite coupling U > Uc, 
since the divergence of Fq is cut-off in this case. In fact. Fig. 2b shows that, for U/t = 4 and 
R = 0.94, the singlet ground state is stable on the 10 x 10 lattice, contrary to the 16 x 16 
cluster. The tendency towards FM grows with increasing cluster size also in the mean-field 
approximation, TMA, and at the Gutzwiller variational level. The lattice-size dependence 
of the difference between the energies of a fully polarized state, and of a closed-shell state 
with minimal spin is shown in Fig. 2c. Our QMC data strongly support the FM state 



for / — s> oo. The qualitative agreement of TMA with the QMC data is striking. TMA 
tends to shghtly overestimate the FM correlations, but the improvement with respect to the 
mean-field approximation and the Gutzwiller wavefunction (see Fig. 1) is evident. 

In order to confirm the prediction of AFM and superconductivity in the remaining part of 
the phase diagram, we have studied the imaginary time dependence of correlation functions 
in the following form: 

For (3 —>■ oo, 0(0) tends to the so called mixed average 0(0) = which is a property of 

both the ground state and the trial wavefunction itself. Instead, for r ~ [3/2, one obtains the 
desired ground-state expectation value of O. Thus, if we choose as a trial wavefunction the 
singlet paramagnetic Gutzwiller wavefunction, i.e. a state without BCS or AFM order, an 
increase of the corresponding correlations is expected by varying r from to (3/2, whenever 
the electron system truly supports some kind of order. The time-dependent functions 0{t) 
defined in (^), directly indicate enhanced or depressed correlations with respect to a reference 
finite-size Gutzwiller state {ipr)- They represent therefore clear fingerprints of the nature of 
the ground state {ipo) even with a short imaginary-time propagation (3, provided the reference 
state \iPt) is sufficiently close to \ipo). A similar idea was used also in Ref. fl^. 



In Fig. 3a we plot the spin-isotropic magnetic structure factor for the AFM wavevector 
Q = (tt, tt) as a function of the lattice size /. Only for / > 12 the AFM correlations are 
clearly enhanced, supporting the existence of a true AFM long range order. Unfortunately, 
the need for /3t > 2 at large / precludes the possibility of a detailed finite-size scaling of 
AFM order. For negative U (Fig. 3b) s-wave superconductivity with an order parameter 
s = I; Z]fc c^,tcLfc_| should take place. The plot of 0(r) where O = s^'s clearly shows this 
tendency even for a short imaginary-time propagation, in agreement with Ref. 0. 

Based on the results of QMC simulations and on the analysis of TMA we conclude 
that, in the t-t' Hubbard model at the VH density, interactions tend to remove the singular 
density of states, via a splitting of the Fermi surface (FM) or the opening of a gap at the 



VH points (AFM and superconductivity). This has been shown numerically at realistic 
coupling strengths U and is expected analytically for all U, if we neglect the change of the 
quasiparticle dispersion. Since we find magnetism at f/ > despite the finite-size cut-offs, 
we expect that it appears, for U large enough, in a finite window of densities around the 
VH filling. Our results for repulsive interactions reduce in the limits R —* and -R ^ 1 
to AFM of the half-filled t' = Hubbard model and to an analogue of the fiat-band FM 
10], respectively. An interesting question arises about the transition from AFM to FM as 
R increases. Since already on the mean field level both FM and AFM susceptibilities are 
suppressed for R ^ 0.55 as compared to their values at R = 0.94 and R = 0.2, respectively, a 
numerical solution of this question will require even bigger lattices than those used here. We 
cannot exclude that close to R ~ 0.55, there is no magnetic ordering and a superconducting 



instability develops via the Kohn-Luttinger effect [15 
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FIGURES 

FIG. 1. The energy per site (in units of t) as a function of m for N = + Ni = 108 electrons. 
Full lines: free electrons (bars), mean-field approximation (triangles), 2°^^ order perturbation 
theory (squares), and T-matrix approximation without self-consistency condition for the electron 
Green's function (circles). Dashed lines: optimal Gutzwiller wavefunction (squares) and QMC data 
at /3t = 2 (triangles). 

FIG. 2. The energy per site (in units of t) (a) of a L = 16 x 16 lattice with N = 108 |^ 
electrons as a function of (3, for A''^ = 55 (total spin 5 = 1, continuous line) and the f3 converged 
one for the FM state (total spin S = 53, dashed line); (b) as a function oira at (3t = 2 for N = 108, 
L = 16 X 16 (upper curve) and N = 44, L = 10 x 10 (lower curve). The optimal g =~ 0.65 (see 
text) is used in all simulations. The error due to the imaginary time discretization (A/3t = 0.2) 
is negligible, as the new estimator ^ considerably reduces both this (~ 10 100 times) and the 
statistical error (3-^4 times) at the optimal g. (c) Difference between the energy per site (in units 
of t) of the minimal-spin [S = 1) state and of the fully polarized state, both at the VH density, as 
a function of the lattice size. 

FIG. 3. QMC simulation for t' /t = 0.1 at the VH density in the singlet sector, (a) Isotropic 
magnetic structure factor S'(7r,7r) = X^ij < Si ■ Sj > g^Qi^i-Rj) at the AF wavevector Q as a 
function of r (see text) for U/t = 4. The curves (guides to the eye) correspond to square lattices 
/ X / with / = 6,8,10,12,16,20 {N„ = 15,27,43,63,115,183) from bottom to top. The optimal 
Gutzwiller parameter g = 0.625. (b) The same plot as in (a) for the s-wave BCS order parameter 
described in the text on a lattice 16 x 16 for U/t = —2. The optimal g = —0.3 and Apt = 0.1. 



10 




R. Hlubina et al. Fig. 1 



0.72 



0.74 



0.76 - 



0.78 



1 \ \ \ \ \ \ \ \ ^ \ \ \ \ ^ \ [ 

FERRO 

SINGLET H 




1^1 A 



0.77 
0.78 
0.79 





1 2 




3 


r(b) u/t=4 ty 

=^^1 — \ — r^h—r 


t-0.47^ 1 

1 — ®^ ~ 
1 1 1 1 - 



0.1 0.2 0.3 0.4 0.5 

m 



R. Hlubina et al. Fig. 2 



o 
Pi 



o 

GO 



0.02 







0.02 



0.04 



1 1 A_i4_-j4— 

- (c) 


^ ' iA ' 


1 1 


- • QMC 






aT-MATRIX 




1 1 1 1 


1 1 1 1 


1 1 







0.05 



0.1 



1/1 



R. Hlubina et al. Fig. 2 




R. Hlubina et al. Fig. 3 



